using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.QuadPack
{
    /// <summary>
    /// DQAGS - Integration over finite intervals.
    /// </summary>
    [Serializable]
    public static class DqagsClass
    {
        /// <summary>
        /// The adaptive strategy compares results of integration over the given interval with the sum of 
        /// results obtained from integration over a bisected interval. Since error estimates are available from 
        /// each regional integration, the region with the largest error is bisected and new results are computed. 
        /// This bisection process is continued until the error is less than the prescribed limit or convergence 
        /// failure is determined.
        /// </summary>
        /// <param name="f">Double precision function to be integrated.</param>
        /// <param name="a">Lower limit of integration.</param>
        /// <param name="b">Upper limit of integration.</param>
        /// <param name="epsabs">The absolute accuracy requested.</param>
        /// <param name="epsrel">The relative accuracy requested.</param>
        /// <param name="abserr">The requested absolute error of the result.</param>
        /// <param name="neval">The needed function evalations.</param>
        /// <param name="ier">The requested errorcode.</param>
        /// <returns>The result of the numerical integration.</returns>
        public static double Dqags(HardRealFunction f, double a, double b, double epsabs, double epsrel,
                                   ref double abserr, ref int neval, ref int ier)
        {
            double abseps = 0;
            double[] alist = new double[QuadPackConstants.LIMIT];
            double area;
            double area1;
            double area12;
            double area2;
            double a1;
            double a2;
            double[] blist = new double[QuadPackConstants.LIMIT];
            double b1;
            double b2;
            double correc = 0;
            double defabs = 0;
            double defab1 = 0;
            double defab2 = 0;
            double dres;
            double[] elist = new double[QuadPackConstants.LIMIT];
            double erlarg = 0;
            double erlast;
            double errbnd;
            double errmax;
            double error1 = 0;
            double error2 = 0;
            double erro12;
            double errsum;
            double ertest = 0;
            double resabs = 0;
            double reseps;
            double result;
            double[] res3la = new double[3];
            double[] rlist = new double[QuadPackConstants.LIMIT];
            double[] rlist2 = new double[52];
            double small = 0;

            int id;
            int ierro;
            int[] iord = new int[QuadPackConstants.LIMIT];
            int iroff1;
            int iroff2;
            int iroff3;
            int jupbnd;
            int k;
            int ksgn;
            int ktmin;
            int last;
            int maxerr;
            int nres;
            int nrmax;
            int numrl2;
            int limit;
            int extrap;
            int noext;

            limit = QuadPackConstants.LIMIT - 1;
            // Test validity of parameters. 	
            ier = 0;
            neval = 0;
            result = 0.0;
            abserr = 0.0;
            alist[0] = a;
            blist[0] = b;
            rlist[0] = 0.0;
            elist[0] = 0.0;
            if ((epsabs < 0.0) && (epsrel < 0.0))
            {
                ier = 6;
            }
            if (ier == 6)
            {
                return result;
            }

            // First approximation to the integral. 
            ierro = 0;
            result = Dqk21Class.G_K21(f, a, b, ref abserr, ref defabs, ref resabs);

            // Test on accuracy. 
            dres = Math.Abs(result);
            errbnd = (((epsabs) > (epsrel * dres)) ? (epsabs) : (epsrel * dres));
            last = 1;
            rlist[0] = result;
            elist[0] = abserr;
            iord[0] = 0;
            if ((abserr <= 100.0 * Double.Epsilon * defabs) && (abserr > errbnd))
            {
                ier = 2;
            }
            if (limit == 0)
            {
                ier = 1;
            }
            if ((ier != 0) || (abserr <= errbnd) && (abserr != resabs) || (abserr == 0.0))
            {
                goto _140;
            }

            // Initialization. 
            rlist2[0] = result;
            errmax = abserr;
            maxerr = 0; // maxerr = 1 
            area = result;
            errsum = abserr;
            abserr = Double.MaxValue;
            nrmax = 0;
            nres = 0; // nres = 0 
            numrl2 = 1; // numrl2 = 2 
            ktmin = 0;
            extrap = QuadPackConstants.FALSE;
            noext = QuadPackConstants.FALSE;
            ierro = 0;
            iroff1 = 0;
            iroff2 = 0;
            iroff3 = 0;
            ksgn = -1;
            if (dres > (1.0 - 50.0 * Double.Epsilon) * defabs)
            {
                ksgn = 1;
            }

            // Main loop. 
            for (last = 1; last <= limit; last++)
            {
                // Bisect the interval with the nrmax-th largest error estimate. 
                a1 = alist[maxerr];
                b1 = 0.5 * (alist[maxerr] + blist[maxerr]);
                a2 = b1;
                b2 = blist[maxerr];
                erlast = errmax;
                area1 = Dqk21Class.G_K21(f, a1, b1, ref error1, ref resabs, ref defab1);
                area2 = Dqk21Class.G_K21(f, a2, b2, ref error2, ref resabs, ref defab2);

                // Improve previous approxminations to integral and error
                //  	and test for accuracy. 
                area12 = area1 + area2;
                erro12 = error1 + error2;
                errsum = errsum + erro12 - errmax;
                area = area + area12 - rlist[maxerr];
                if ((defab1 == error1) || (defab2 == error2))
                {
                    goto _15;
                }
                if ((Math.Abs(rlist[maxerr] - area12) > 1.0e-5 * Math.Abs(area12)) || (erro12 < .99 * errmax))
                {
                    goto _10;
                }
                if (extrap != 0)
                {
                    iroff2++;
                }
                else
                {
                    iroff1++;
                }
                _10:
                if ((last > 9) && (erro12 > errmax)) // last > 10 
                {
                    iroff3++;
                }
                _15:
                rlist[maxerr] = area1;
                rlist[last] = area2;
                errbnd = (((epsabs) > (epsrel * Math.Abs(area))) ? (epsabs) : (epsrel * Math.Abs(area)));

                // Test for roundoff error and eventually set error flag. 
                if (((iroff1 + iroff2) >= 10) || (iroff3 >= 20))
                {
                    ier = 2;
                }
                if (iroff2 > 5)
                {
                    ier = 3;
                }

                // Set error flag in the case that the number of subintervals
                //	equals limit. 
                if (last == limit) // last == limit 
                {
                    ier = 1;
                }

                // Set error flag in the case of bad integrand behavior at some
                //	points in the integration range. 
                if ((((Math.Abs(a1)) > (Math.Abs(b2))) ? (Math.Abs(a1)) : (Math.Abs(b2))) <=
                    (1.0 + 1000.0 * Double.Epsilon) * (Math.Abs(a2) + 1000.0 * Double.MinValue))
                {
                    ier = 4;
                }

                // Append the newly-created intervals to the list. 
                if (error2 > error1)
                {
                    goto _20;
                }
                alist[last] = a2;
                blist[maxerr] = b1;
                blist[last] = b2;
                elist[maxerr] = error1;
                elist[last] = error2;
                goto _30;
                _20:
                alist[maxerr] = a2;
                alist[last] = a1;
                blist[last] = b1;
                rlist[maxerr] = area2;
                rlist[last] = area1;
                elist[maxerr] = error2;
                elist[last] = error1;

                // Call dqsort to maintain the descending ordering in the list of error
                //	estimates and select the subinterval with nrmax-th largest
                //	error estimate (to be bisected next). 
                _30:
                DqsortClass.Dqsort(limit, last, ref maxerr, ref errmax, elist, iord, ref nrmax);
                if (errsum <= errbnd)
                {
                    goto _115;
                }
                if (ier != 0)
                {
                    goto _100;
                }
                if (last == 1) // last == 2 
                {
                    goto _80;
                }
                if (noext != 0) // goto 90 
                {
                    goto _90;
                }
                erlarg -= erlast;
                if (Math.Abs(b1 - a1) > small)
                {
                    erlarg += erro12;
                }
                if (extrap != 0)
                {
                    goto _40;
                }

                // Test whether the interval to be bisected next is the smallest interval. 
                if ((Math.Abs(blist[maxerr] - alist[maxerr])) > small)
                {
                    goto _90; // goto 90 
                }
                extrap = QuadPackConstants.TRUE;
                nrmax = 1; // nrmax = 2 
                _40:
                if ((ierro == 3) || (erlarg <= ertest))
                {
                    goto _60;
                }

                // The smallest interval has the largest error. Before bisecting, decrease
                //	the sum of the erors over the larger intervals (erlarg) and
                //		perform extrapolation.) 
                id = nrmax;
                jupbnd = last;
                if (last > (2 + limit / 2))
                {
                    jupbnd = limit + 3 - last;
                }
                for (k = id; k <= jupbnd; k++)
                {
                    maxerr = iord[nrmax];
                    errmax = elist[maxerr];
                    if (Math.Abs(blist[maxerr] - alist[maxerr]) > small)
                    {
                        goto _90; // goto 90 
                    }
                    nrmax++;
                }

                // Perform extrapolation. 
                _60:
                numrl2++;
                rlist2[numrl2] = area;
                reseps = DqextClass.Dqext(ref numrl2, rlist2, ref abseps, res3la, ref nres);
                ktmin++;
                if ((ktmin > 5) && (abserr < 1.0e-3 * errsum))
                {
                    ier = 5;
                }
                if (abseps >= abserr)
                {
                    goto _70;
                }
                ktmin = 0;
                abserr = abseps;
                result = reseps;
                correc = erlarg;
                ertest = (((epsabs) > (epsrel * Math.Abs(reseps))) ? (epsabs) : (epsrel * Math.Abs(reseps)));
                if (abserr <= ertest)
                {
                    goto _100;
                }

                // Prepare bisection of the smallest interval. 
                _70:
                if (numrl2 == 0)
                {
                    noext = QuadPackConstants.TRUE;
                }
                if (ier == 5)
                {
                    goto _100;
                }
                maxerr = iord[0];
                errmax = elist[maxerr];
                nrmax = 0;
                extrap = QuadPackConstants.FALSE;
                small = small * 0.5;
                erlarg = errsum;
                goto _90; // goto 90 
                _80:
                small = Math.Abs(b - a) * 0.375;
                erlarg = errsum;
                ertest = errbnd;
                rlist2[1] = area;
                _90:
                ;
            } // 90: 
            _100:
            if (abserr == Double.MaxValue)
            {
                goto _115;
            }
            if ((ier + ierro) == 0)
            {
                goto _110;
            }
            if (ierro == 3)
            {
                abserr += correc;
            }
            if (ier == 0)
            {
                ier = 3;
            }
            if ((result != 0.0) && (area != 0.0))
            {
                goto _105;
            }
            if (abserr > errsum)
            {
                goto _115;
            }
            if (area == 0.0)
            {
                goto _130;
            }
            goto _110;
            _105:
            if (abserr / Math.Abs(result) > errsum / Math.Abs(area))
            {
                goto _115;
            }

            // Test on divergence. 
            _110:
            if ((ksgn == -1) &&
                ((((Math.Abs(result)) > (Math.Abs(area))) ? (Math.Abs(result)) : (Math.Abs(area))) <= defabs * .01))
            {
                goto _130;
            }
            if ((0.01 > result / area) || (result / area > 100.0) || (errsum > Math.Abs(area)))
            {
                ier = 6;
            }
            goto _130;

            // Compute global integral. 
            _115:
            result = 0.0;
            for (k = 0; k <= last; k++)
            {
                result += rlist[k];
            }
            abserr = errsum;
            _130:
            if (ier > 2)
            {
                (ier)--;
            }
            _140:
            neval = 42 * last - 21;
            return result;
        }
    }
}